/*
 *  flexfem.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 02/04/09.
 * 
 *  Modified by Markus Lutz on 08/10/09
 */

#include "flexfem.h"

FlexFem::FlexFem(ifstream &in, string cname, MecSystem *pSys)
:	Body(in,cname,pSys),
	meshFileName(),
	nNodes(),
	nElements(),
//	vPointNames(),
//	connect(),
	connectSurface(), // edited
	vCoord(),
	vertexTag(),
	nVertex(),
	matName("")
{	
	dim = parentSystem->dim;	/** Must be set always.  dim belongs to the mother class Body */
	type = "flexfem";
	
	vCoord.resize(dim);
	
 	ifstream ini;
	char keyword[20];
	string auxs;				/** Aux string */
	vector<vector<string> > connect;
	char contactFileName[30];
	//int numContact;
	string line;
	vector<string> vPNames;
	int i;
	//ifstream inContact;

	while (in >> keyword && strcmp(keyword,"ENDBODY") )
	{
		if (!strcmp(keyword,"MESHFILE"))
		{
			/** keyword has MESHFILE */
			in >> meshFileName;	/** Reads the name of the file */
		
			in >> auxs;				/** Reads the MAT label */
			in >> matName;			/** Reads the material name */
			ini.open(meshFileName);
			// gmsh format
			readGmsh(ini,connect);
			getBoundary(connect); // edit
			ini.close();
		
			// NOT fully implemented yet
			// createSurface();
		}
		if (!strcmp(keyword,"CONTACTPOINTS"))
		{
			in >> auxs;
			if (!auxs.compare("all"))	contactAllOption = true;
			if (!auxs.compare("file"))
			{
				in >> contactFileName;
				contactFileOption = true;
			}
			if (!auxs.compare("names"))
			{
				contactNamesOption = true;
				getline(in,line);
				istringstream iss(line);
				while (iss >> keyword )
				{
					auxs = keyword;
					vPNames.push_back(auxs);
					//cout << "keyword is " << keyword << endl;
				}				
			}
		}				
		if(parentSystem->getProcNum()==parentSystem->getRankNum()) cout << "	Body " << name << " succesfully read \n";
	}

	// If nothing is said about contact, take by default the "all" option
	if (contactAllOption==false && contactFileOption==false && contactNamesOption==false)	contactAllOption = true;
	
	//cout << "contactAllOption = " << contactAllOption << endl;
	//cout << "contactFileOption = " << contactFileOption << endl;
	//cout << "contactNamesOption = " << contactNamesOption << endl;
	if (contactAllOption==true)
	{  					
		vPointNamesContact = vPointNamesSurface;
	}
	if (contactFileOption==true)
	{
		readContactFromFile(contactFileName);
		/*
		inContact.open(contactFileName);
		inContact >> numContact;
		for (i=0; i<numContact; ++i)
		{
			inContact >> auxs;
			vPointNamesContact.push_back(auxs);
		}
		inContact.close();
		*/
	}
	if (contactNamesOption==true)
	{
		for (i=0; i<vPNames.size(); ++i)	vPointNamesContact.push_back(vPNames[i]);
	}
	//cout << "vPointNamesContact " << endl;
	//printSTLVector(vPointNamesContact);
}

// ***

void FlexFem::setPTPointsComp()
{
	//cout << "EN bodyyyyyy " << name << endl;
	int j, k;
	vector<string> auxs;
	map<string, Point*, less<string> > maux;
	pair<string,Point*> mauxEntry;
	pair<string,int> sharedEntry;
	
	for(j=0; j<vComponents.size(); j++)		
	//Loop on the components
	{
		auxs = (vComponents[j]->getPointNames());
		for(k=0; k<auxs.size(); k++)		
		// Loop on the names of the points within each component
		{
//			cout << "puntoooo " << auxs[k] << endl;
			// First, check if the point has been renamed
//			if( (parentSystem->checkIfPointChanged(auxs[k]))==true)
//			{
//				cout << "He encontrado un punto cambiado, es el " << auxs[k] << endl;
//			}
			/** If the point it is not at the local process ... */
			if( (parentSystem->checkIfPointProcess(auxs[k]))==false )
			{
				 //cout << "Punto " << auxs[k] << " es remoto, estoy en el sistema " << parentSystem->name << endl;
				/** ...  create it at root system, copying the original point from the remote process */
				parentSystem->createPointAtRoot( parentSystem->getPTPoint(auxs[k]) );
				/** ...  and add the name of the shared point at the remote process */
				parentSystem->addRemoteSharedPoint(auxs[k],parentSystem->getProcNum());
			}
			mauxEntry.first = auxs[k];
			mauxEntry.second = (parentSystem->getPTPointLocal(auxs[k]));
			maux.insert(mauxEntry);
		}
		vComponents[j]->setMPoints(maux);
		maux.clear();
	}
}

// ***

void FlexFem::setPTPointsLoad()
{
	int j, k;
	vector<string> auxs;
	map<string, Point*, less<string> > maux;
	pair<string,Point*> mauxEntry;
	pair<string,int> sharedEntry;
	
	for(j=0; j<vLoads.size(); j++)		
	//Loop on the Loads
	{
		auxs = (vLoads[j]->getPointNames());
		for(k=0; k<auxs.size(); k++)		
		// Loop on the names of the points within each component
		{
			/** If the point it is not at the local process ... */
			if( (parentSystem->checkIfPointProcess(auxs[k]))==false )
			{
				/** ...  create it at root system, copying the original point from the remote process */
				parentSystem->createPointAtRoot( parentSystem->getPTPoint(auxs[k]) );
				/** ...  and add the name of the shared point at the remote process */
				parentSystem->addRemoteSharedPoint(auxs[k],parentSystem->getProcNum());
			}
			mauxEntry.first = auxs[k];
			mauxEntry.second = (parentSystem->getPTPointLocal(auxs[k]));
			maux.insert(mauxEntry);
		}
		vLoads[j]->setMPoints(maux);
		maux.clear();
	}
}

// ***

void FlexFem::setPPointSurface()
{

	// Build the vector of pointers to the surface points
	vPPointSurface.resize(vPointNamesSurface.size());
	
	int i;
	string auxs;	
	for(i=0;i<vPointNamesSurface.size();i++)
	{
		auxs = vPointNamesSurface[i];
		addSystemName(auxs,parentSystem->name);
		vPPointSurface[i] = parentSystem->getPTPointLocal( auxs ); 
	}
}

// ***

void FlexFem::setMatAtComp()
{
	int i;
	string matName;
	Material* mat;
	for(i=0; i<vComponents.size(); i++)
	{
		matName = vComponents[i]->getMatName();
		mat = parentSystem->getPTMaterial(matName);		
		vComponents[i]->setMatAtComp(mat);
	}
}

// ***

void FlexFem::getBoundary(vector<vector<string> >& connect)
{
	int i;
	
// edit1: starts here ....

// this works for cubes only, further geometry has to added...

//*********************************************
// converting the cube into surfacees
//*********************************************

	switch (vertexTag) // Element tag (2=triangle, 3=quadrangle, 4=tetrahedron, 5=hexahedron)
	{
		case 2: {break;}
		
		case 3: {break;} 

		case 4: {break;}

		case 5: // hexadron
		{
		
			vector <string> surface(5,"0");		// tempory vector for the cuadro + a identifier 

			// now the cube is reconstructed by 6 surfaces
			for(i = 0; i < nElements;  i += 1)	// loop over all Elements
			{	
		
				// first surface
				surface[0] = connect[i][0];
				surface[1] = connect[i][1];
				surface[2] = connect[i][2];
				surface[3] = connect[i][3];
				surface[4] = connect[i][0]+connect[i][1]+connect[i][2]+connect[i][3]; // identifier

				connectSurface.push_back(surface);


				// second surface
				surface[0] = connect[i][4];
				surface[1] = connect[i][0];
				surface[2] = connect[i][3];
				surface[3] = connect[i][7];
				surface[4] = connect[i][4]+connect[i][0]+connect[i][3]+connect[i][7]; // identifier

				connectSurface.push_back(surface);


				// third surface
				surface[0] = connect[i][4];
				surface[1] = connect[i][5];
				surface[2] = connect[i][1];
				surface[3] = connect[i][0];
				surface[4] = connect[i][4]+connect[i][5]+connect[i][1]+connect[i][0]; // identifier

				connectSurface.push_back(surface);


				// fourth surface
				surface[0] = connect[i][4];
				surface[1] = connect[i][5];
				surface[2] = connect[i][6];
				surface[3] = connect[i][7];
				surface[4] = connect[i][4]+connect[i][5]+connect[i][6]+connect[i][7]; // identifier

				connectSurface.push_back(surface);


				// fifth surface
				surface[0] = connect[i][5];
				surface[1] = connect[i][1];
				surface[2] = connect[i][2];
				surface[3] = connect[i][6];
				surface[4] = connect[i][5]+connect[i][1]+connect[i][2]+connect[i][6]; // identifier

				connectSurface.push_back(surface);


				// sixth surface
				surface[0] = connect[i][7];
				surface[1] = connect[i][6];
				surface[2] = connect[i][2];
				surface[3] = connect[i][3];
				surface[4] = connect[i][7]+connect[i][6]+connect[i][2]+connect[i][3]; // identifier

				connectSurface.push_back(surface);
		
			}
			break;
		}
		default:	{cout << "vertexTag" << vertexTag << "unknown\n"; break;}
	}	


//*******************************************************
// reductant surfaces are deleted
//*******************************************************

	// vectors	
	vector< vector<string> > save_vector = connectSurface; 	// vector which saves the original structure of the surfacees
	vector< vector<string> > temp_vector;			// vector for the elements of the surface
	vector<string> match;					// vector for the element identifiers of the surface


	// the nodes in all elements are sorted, last element of the vector is not considered
	for (i = 0; i<connectSurface.size();i++){
		sort (connectSurface[i].begin(), connectSurface[i].end()-1);	
	}

	
	// now the vectors containing the nodes are sorted
	sort(connectSurface.begin(),connectSurface.end());


	// as the vectors are sorted now, only the neighbour of each vector
	// has to be compared for equality

	
  	for (i = 0; i<connectSurface.size()-1;i++)	// loop over all vectors
	{
		if (!equal( connectSurface[i].begin(), connectSurface[i].end()-1, connectSurface[i+1].begin())){
			match.push_back(connectSurface[i][4]);
		}
		else {i += 1;}
	}

	// the last surface is a special case...
	if (!equal( connectSurface[connectSurface.size()-1].begin(), connectSurface[connectSurface.size()-1].end()-1, connectSurface[connectSurface.size()-2].begin()))
	{
		match.push_back(connectSurface[connectSurface.size()-1][4]);
	}


	sort(match.begin(),match.end());	// arrange vector

	for (i = 0; i< connectSurface.size(); i++ ) // loop over all surfaces
	{	
		if (binary_search(match.begin(),match.end(),save_vector[i][4])){	// elements of the surface are compared to the vector containing all
			temp_vector.push_back(save_vector[i]);				// "new" elements, if found
		}
	}
	connectSurface.resize(temp_vector.size());
	connectSurface = temp_vector;		// vector with new elements of the surface

//*******************************************************
// reductant nodes are deleted
//*******************************************************
		
	for (i = 0; i<connectSurface.size();i++) 			// make big vector with all nodes
	{
		vPointNamesSurface.push_back(connectSurface[i][0]);
		vPointNamesSurface.push_back(connectSurface[i][1]);
		vPointNamesSurface.push_back(connectSurface[i][2]);
		vPointNamesSurface.push_back(connectSurface[i][3]);
	}
		
	sort(vPointNamesSurface.begin(),vPointNamesSurface.end()); 	// arrange nodes 	
		
	vector <string>::iterator it;
	it = unique (vPointNamesSurface.begin(), vPointNamesSurface.end()); 	// delete all multiple nodes in vector
  	vPointNamesSurface.resize( it - vPointNamesSurface.begin() ); 		// resize node vector

// and edit1 ends here
}

// ***

void FlexFem::readGmsh(ifstream &in, vector<vector<string> >& connect)
{
	string auxs;
	vector<string> vVertex;		/** Auxiliary vector with the names of the vertexes in one element */
	int auxi1, auxi2, auxi3;

	std::string version_number, file_type, data_size;
	int number_of_tags, tag;

	int i;
	int j;
	double coo;
//	vector<string> auxsv;		/** Aux string vector. Contains names of vertexes. */
//	auxsv.resize(8);
	vector<string> vPointNames;

	// $MeshFormat
	in >> auxs;
	in >> version_number >> file_type >> data_size;
	in >> auxs;
	// $EndMeshFormat

	in >> auxs;		/** $Nodes */
	in >> nNodes;
	vPointNames.resize(nNodes);

	for (i=0; i<nNodes; ++i)	// Read coordinates
	{
		in >> auxs;
		vPointNames[i] = auxs;				// Name of the node
		addSystemName(auxs,parentSystem->name);			/** System name is added to the point name */

		in >> coo;			// First coordinate
		vCoord[0] = coo;
		in >> coo;			// Second coordinate
		vCoord[1] = coo;
		in >> coo;
		if (dim==3) vCoord[2] = coo;	// Third coordinate.
		
		// It is assumed that nodes do not exist previously.  Their names must be added at parentSystem
		parentSystem->addPointName(auxs,vCoord);
	}
	in >> auxs;		/** $EndNodes */
	in >> auxs;		/** $Elements */
	in >> nElements;
	connect.resize(nElements);

	for (i=0; i<nElements; ++i)
	{
		in >> auxi1;	// Element number (not used)
		in >> vertexTag;// Element tag

		if (i==0)
		{
			nVertex = 0;

			// 2 = triangle
			if (vertexTag == 2)
				nVertex = 3;

			// 3 = quadrangle
			if(vertexTag == 3)
				nVertex = 4;

			// 4 = tetrahedron
			if(vertexTag == 4)
				nVertex = 4;

			// 5 = hexahedron
			if(vertexTag == 5)
				nVertex = 8;

			//Element Type not used in Eppi.
			if(nVertex == 0)
				cout << "**ERROR: Mesh non-consistent with dimension\n";

			vVertex.resize(nVertex);
//			auxsv.resize(nVertex);
			for (j = 0; j < nElements; ++j)
			{
				connect[j].resize(nVertex);
			}
		}

		in >> number_of_tags;
		//cout << number_of_tags;
		for(j = 0; j < number_of_tags; ++j)
		{
			in >> tag;
		}

		for (j=0; j<8; ++j)	// Maximum number of vertexes is 8
		{
			in >> auxs;	// First vertex
			connect[i][j] = auxs;
			addSystemName(auxs,parentSystem->name);	/** System name is added to the point name */
			vVertex[j] = auxs;
			if (j==2 && nVertex<4) break;	// Triangles
			if (j==3 && nVertex==4) break;	// Quadrilaterals or tetrahedra
		}
		if (vertexTag==5) vComponents.push_back(new Brick(vVertex,matName,dim));	// Only hexahedron. Other elements TO BE DONE
	}
	in >> auxs;		/** $EndElements */
}

void FlexFem::writeMSH(ofstream &fGraph)
{
	int i;
	int j;
	int nt;
	
	string numStep = "0"; 	// edit: current timestep 
	
	fGraph << "$MeshFormat\n";
	fGraph << "2 0 8\n";
	fGraph << "$EndMeshFormat\n";
	fGraph << "$Nodes\n";
	fGraph << vPPointSurface.size() <<  "\n"; // edit
	
	for (i=0; i<vPPointSurface.size(); ++i)
	{
		nt = vPPointSurface[i]->globalTNum.size();	// Number of translational dofs
		fGraph << vPointNamesSurface[i]; // edit 
		for (j=0; j<nt; ++j)
		{
			fGraph << "\t" << parentSystem->getDof( vPPointSurface[i]->globalTNum[j] ); 
		}
		if (nt<3) fGraph << "\t" << 0;	// Third coordinate, just in case is 2D
		fGraph << "\n";
	}
	fGraph << "$EndNodes\n";
	
	
	fGraph << "$Elements\n";
	fGraph << connectSurface.size() << "\n";
	for (i=0; i<connectSurface.size(); ++i)
	{
		fGraph << i << "\t" << "3 " << "\t" << "3 " << "7 " << "6 "<< "0 ";  // edit
		// fGraph << i << "\t" << vertexTag << "\t" << "3 " << "7 " << "6 "<< "0 ";	
		for (j=0; j<4; ++j) // edit
			// for (j=0; j<nVertex; ++j)
		{
			fGraph << "\t" << connectSurface[i][j]; // edit
			// fGraph << "\t" << connect[i][j];
		}
		fGraph << "\n";
	}
	fGraph << "$EndElements\n";
	
	// get current time step from the outputfile name stored for that process
	// the graphic file name has to be: xxxxx-number.dat, number will be extracted from the file name, where "-" and "." are only used one time
	
	string outGraphFile(parentSystem->getOutGraphFile());			// get name of the graphic file
	outGraphFile = outGraphFile.substr(0,outGraphFile.find_last_of(".")); 	// get first part of the filename until ".(dat)"
	numStep = outGraphFile.substr(outGraphFile.find_last_of("-")+1);	// isolate current step
	
	
	fGraph << "$NodeData\n"; 
	
	fGraph << "1\n"; 			// one string follows
	fGraph << "\"Title or new name\"\n"; 	// name of the "view"
	
	fGraph << "1\n";			// one real tag follows
	fGraph << numStep << "\n";		// time value // can be different, should be changed to real time value !!!!
	
	fGraph << "3\n";			// three integer will follow
	fGraph << numStep << "\n";		// time step
	fGraph << "1\n";				// type: 1 = scalar field,...			 
	fGraph << vPointNamesSurface.size() <<"\n";	// number of associated nodal values
	
	// now start the nodal values
	for (int i = 0; i < vPointNamesSurface.size(); i++)
	{
		fGraph << vPointNamesSurface[i] << "\t" << "0.1" << "\n"; 	// "0.1" should be replaced by displacement or tension
	}
	
	fGraph << "$EndNodeData\n"; 
	
	
	/*****************************************************/
	
	/** file with the gmsh command for the post processing animation*/ 
	// this is only a test a should be deleted later
	
	ofstream gmshcom;
	
	if (!(numStep.compare("0")))
	{
		gmshcom.open("gmsh_com.txt");
		gmshcom.close();
	}
	//double  test;
	//cin >> test;
	
	gmshcom.open("gmsh_com.txt",ios::app);
	
	if (!(numStep.compare("0")))
	{
		gmshcom << "gmsh -combine ";
	}
	
//	string outGraphFilex(parentSystem->getOutGraphFile());	
	outGraphFilex = parentSystem->getOutGraphFile();	
	gmshcom << outGraphFilex << " -new ";
	
	
	gmshcom.close()	;	
}

// ***

void FlexFem::writeVTU(ofstream &fGraph, bool &plotDispFlag, bool &plotProcFlag, bool& plotStressFlag)
{
	int i;
	int j;
	int nt;
	
/*
	string s1, s2, s3;

	outGraphFilex = parentSystem->getOutGraphFile();
	s1 = outGraphFilex.substr(0,outGraphFilex.find_last_of("."));	// Remove the extension
	s2 = s1.substr(5);	// Remove Graph (it has 5 letters)
	s3 = name + s2 + ".vtu";
*/	
//	fGraph.open(outGraphFile,ios::app);

	localFGraph << "<?xml version=\"1.0\"?>\n";
	localFGraph << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
	localFGraph << "	<UnstructuredGrid>\n";
	localFGraph << "		<Piece NumberOfPoints=\"" << vPPointSurface.size() << "\" NumberOfCells=\"" << connectSurface.size() << "\">\n";
	localFGraph << "			<Points>\n";
	localFGraph << "				<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">";
	for(i=0; i<vPPointSurface.size(); ++i)
	{
		localFGraph << "\n";
		localFGraph << "				";
		nt = vPPointSurface[i]->globalTNum.size();	// Number of translational dofs
		for (j=0; j<nt; ++j)
		{
			localFGraph << parentSystem->getDof( vPPointSurface[i]->globalTNum[j] ) << "\t"; 
		}
		if (nt<3) localFGraph << 0.0;	// Third coordinate, just in case is 2D
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph << "			</Points>\n";
	localFGraph<< "			<Cells>\n";
	localFGraph << "				<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">";
	
	// WARNING: Note that the last entry of each row of connectSurface is an internal identifier, and should not be read
	int css0 = connectSurface[0].size()-1;	
	for (i=0; i<connectSurface.size(); ++i)
	{
		localFGraph << "\n";
		localFGraph << "					";
		for (j=0; j<css0; ++j) // edit
		{
			localFGraph << pNamePNumber[ connectSurface[i][j] ] << "\t";
		}
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph << "				<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">";
	localFGraph << "\n";
	localFGraph << "					";
	for (i=0; i<connectSurface.size(); ++i)
	{
		localFGraph << css0*(1+i) << "\t";
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph << "				<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">";
	localFGraph << "\n";
	localFGraph << "					";
	for (i=0; i<connectSurface.size(); ++i)
	{
		localFGraph << 9 << "\t";
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph<< "			</Cells>\n";

	if (plotDispFlag==true || plotProcFlag==true || plotStressFlag==true)	
	localFGraph << "			<PointData Scalars=\"Scalars\">\n";

	// Displacements
	if (plotDispFlag==true)
	{
		double disp = 0.0;
		// Components
		for (j=0; j<dim; ++j)
		{
			localFGraph << "				<DataArray type=\"Float64\" Name=\"Displacement-" << j+1 << "\" format=\"ascii\">";
			localFGraph << "\n";
		
			for (i=0; i<vPPointSurface.size(); ++i)
			{
				disp = parentSystem->getDof( vPPointSurface[i]->globalTNum[j] );
				disp -= vPPointSurface[i]->r[j];
				localFGraph << disp << "\t"; 
				//localFGraph << 0.4 << "\n";
			}
			localFGraph << "\n";
			localFGraph << "				</DataArray>\n";
		}
		// Magnitude
		localFGraph << "				<DataArray type=\"Float64\" Name=\"Displacement-T\" format=\"ascii\">";
		localFGraph << "\n";
		for (i=0; i<vPPointSurface.size(); ++i)
		{
			disp = 0.0;
			for (j=0; j<dim; ++j)
			{
				disp +=  pow(parentSystem->getDof( vPPointSurface[i]->globalTNum[j] ) - vPPointSurface[i]->r[j],2);
			}
			localFGraph << sqrt(disp) << "\t";
		}
		localFGraph << "\n";
		localFGraph << "				</DataArray>\n";
	}
	// Process
	if (plotProcFlag==true)
	{
		localFGraph << "				<DataArray type=\"Float64\" Name=\"Process\" format=\"ascii\">";
		localFGraph << "\n";
		for (i=0; i<vPPointSurface.size(); ++i)
		{
			localFGraph << parentSystem->getProcNum() << "\t";
		}
		localFGraph << "\n";
		localFGraph << "				</DataArray>\n";
	}
	// Stress
	if (plotStressFlag==true)
	{
		// Components
		for (j=0; j<dim; ++j)
		{
			localFGraph << "				<DataArray type=\"Float64\" Name=\"Stress-" << j+1 << "\" format=\"ascii\">";
			localFGraph << "\n";
			
			for (i=0; i<vPPointSurface.size(); ++i)
			{
				localFGraph << parentSystem->getStressValue( 0, vPPointSurface[i]->globalTNum[j] ) << "\t"; 
			}
			localFGraph << "\n";
			localFGraph << "				</DataArray>\n";
		}
		localFGraph << "				<DataArray type=\"Float64\" Name=\"Stress-" << 12 << "\" format=\"ascii\">";
		localFGraph << "\n";
		
		for (i=0; i<vPPointSurface.size(); ++i)
		{
			localFGraph << parentSystem->getStressValue( 1, vPPointSurface[i]->globalTNum[0] ) << "\t"; 
		}
		localFGraph << "\n";
		localFGraph << "				</DataArray>\n";
		localFGraph << "				<DataArray type=\"Float64\" Name=\"Stress-" << 13 << "\" format=\"ascii\">";
		localFGraph << "\n";
		
		for (i=0; i<vPPointSurface.size(); ++i)
		{
			localFGraph << parentSystem->getStressValue( 1, vPPointSurface[i]->globalTNum[1] ) << "\t"; 
		}
		localFGraph << "\n";
		localFGraph << "				</DataArray>\n";
		localFGraph << "				<DataArray type=\"Float64\" Name=\"Stress-" << 23 << "\" format=\"ascii\">";
		localFGraph << "\n";
		
		for (i=0; i<vPPointSurface.size(); ++i)
		{
			localFGraph << parentSystem->getStressValue( 1, vPPointSurface[i]->globalTNum[2] ) << "\t"; 
		}
		localFGraph << "\n";
		localFGraph << "				</DataArray>\n";
		// Principal stresses
		for (j=0; j<dim; ++j)
		{
			localFGraph << "				<DataArray type=\"Float64\" Name=\"PStress-" << j+1 << "\" format=\"ascii\">";
			localFGraph << "\n";
			
			for (i=0; i<vPPointSurface.size(); ++i)
			{
				localFGraph << parentSystem->getStressValue( 2, vPPointSurface[i]->globalTNum[j] ) << "\t"; 
			}
			localFGraph << "\n";
			localFGraph << "				</DataArray>\n";
		}
		// Von Mises
		localFGraph << "				<DataArray type=\"Float64\" Name=\"Von Mises\" format=\"ascii\">";
		localFGraph << "\n";
		
		for (i=0; i<vPPointSurface.size(); ++i)
		{
			localFGraph << parentSystem->getStressValue( 3, vPPointSurface[i]->globalTNum[0] ) << "\t"; 
		}
		localFGraph << "\n";
		localFGraph << "				</DataArray>\n";
	}
	if (plotDispFlag==true || plotProcFlag==true || plotStressFlag==true)
	localFGraph << "			</PointData>\n";
	
	localFGraph << "		</Piece>\n";
	localFGraph << "	</UnstructuredGrid>\n";
	localFGraph << "</VTKFile>\n";	
}

// ***

void FlexFem::writeOOGL(ofstream &fGraph)
{}

// ***

void FlexFem::createSurface()
{
	vector<vector<double> > vlp;	// local points
	vector<vector<int> > vle;		// local elements
	int i;
	int j;
	int nt;
	vector<double> auxd;
	
	for(i=0; i<vPPointSurface.size(); ++i)
	{
		nt = vPPointSurface[i]->globalTNum.size();	// Number of translational dofs
		auxd.resize(nt);
		for (j=0; j<nt; ++j)
		{
			auxd.push_back( parentSystem->getDof( vPPointSurface[i]->globalTNum[j] ) ); 
		}
		vlp.push_back(auxd);
		auxd.clear();
	}
	
	vector<int> vElement;
	
	// WARNING: Note that the last entry of each row of connectSurface is an internal identifier, and should not be read
	int css0 = connectSurface[0].size()-1;	
	vElement.resize(css0);
	for (i=0; i<connectSurface.size(); ++i)
	{
		for (j=0; j<css0; ++j)
		{
			vElement.push_back( pNamePNumber[ connectSurface[i][j] ] );
		}
		vle.push_back(vElement);
		vElement.clear();
	}
	string kk = "kk";
	surface = new SurfaceMesh(kk,vlp,vle);
}

// ***

void FlexFem::readContactFromFile(char* cfn)
{
	int i;
//	int numContact;
//	string auxs;
	ifstream inContact;
	string line;
	vector<string> linesFile;
	vector<double> vauxd(3);
	double auxd;
	string nodeName;	
	
	inContact.open(cfn);
	
	while (inContact.good())
	{
		getline(inContact, line);
		linesFile.push_back(line);
	}
	
	//cout << "linesFile is " << linesFile.size() <<  endl;
	//printSTLVector(linesFile);
	//cout << linesFile[0] << endl;
	
	if (!linesFile[0].compare("NODES"))
	{
		for (int i = 1; i < linesFile.size()-1; ++i)
		{
			istringstream iss(linesFile[i], istringstream::in);
			
			iss >> nodeName;
			vPointNamesContact.push_back(nodeName);
		}
	}
	/*
	if (!linesFile[0].compare("COORD"))
	{
		for (int i = 1; i < linesFile.size(); ++i)
		{
			istringstream iss(linesFile[i], istringstream::in);
			
			for (int j = 0; j < 3; ++j)
			{
				iss >> vauxd[j];
			}
			//vContactPoints.push_back(vauxd);
		}			
	}
	*/
	inContact.close();
}

// ***

void FlexFem::writeFileTemp(const double& time, ofstream & outF)
{}

// ***

void FlexFem::setIniPosBody(vector<double>& iniValues)
{}
